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We introduce and study numerically a directed two-dimensional sandpile automaton with probabilis- 
tic toppling (probability parameter p) which provides a good laboratory to study both self-organized 
criticality and the far-from-equilibrium phase transition. In the limit p = 1 our model reduces to the 
critical height model in which the self-organized critical behavior was found by exact solution [D. 
Dhar and R. Ramaswamy, Phys. Rev. Lett. 63, 1659 (1989)]. For < p < 1 metastable columns 
of sand may be formed, which are relaxed when one of the local slopes exceeds a critical value a c . 
By varying the probability of toppling p we find that a continuous phase transition occurs at the 
critical probability p c , at which the steady states with zero average slope (above p c ) are replaced by 
states characterized by a finite average slope (below p c ). We study this phase transition in detail by 
introducing an appropriate order parameter and the order-parameter susceptibility \- in a certain 
range of p < 1 we find the self-organized critical behavior which is characterized by nonuniversal 
p— dependent scaling exponents for the probability distributions of size and length of avalanches. 
We also calculate the anisotropy exponent £ and the fractal dimension df of relaxation clusters in 
the entire range of values of the toppling parameter p. We show that the relaxation clusters in our 
model are anisotropic and can be described as fractals for values of p above the transition point. 
Below the transition they are isotropic and compact. 
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I. INTRODUCTION 



model displays SOC for certain values of p < 1. 



The term " self-organized criticality" (SOC) (!]] refers to 
certain driven dissipative systems which organize them- 
selves into a steady state far away from equilibrium. The 
critical steady state appears as an attractor of the dy- 
namics, and thus no fine tuning of an external parame- 
ter is necessary to achieve criticality. The major aspect 
of the dynamics of SOC systems is the "separation of 
time scales", i.e. the system relaxes infinitely fast com- 
pared to the time between two successive perturbations. 
From this point of view the relaxation processes may 
be considered as consisting of a series of discrete events 
(avalanches) and one investigates the probability distri- 
bution of their spatial properties, for instance the size of 
an avalanche. The self-organized critical state is char- 
acterized by a power-law behavior of these probability 
distributions. For detailed discussions of the definition 
of SOC we refer to gi. 

Sandpile automata are well known prototype 

models exhibiting SOC. By adding a particle from the 
outside, a sandpile is perturbed and the perturbation 
may lead to instabilities at neighboring sites if the local 
height of sand exceeds a critical value h c (critical height 
model (CHM) 0Q), or if the local slopes exceed a critical 
value a c (critical slope model (CSM) ||-0]). The steady 
state of the system in the case of the CHM is accompa- 
nied with a state of a zero average slope. A finite average 
slope characterizes a CSM. 

In Section II we introduce a two dimensional directed 
sandpile model with stochastic dynamics, in which the 
updating rules may be tuned continuously by varying 
a parameter p, representing probability of toppling, be- 
tween CHM (for p = 1) and CSM (for p = 0) §. Our 
numerical results for p < 1 show that a noncquilibrium 
phase transition takes place at the critical value p c , which 
is characterized by a continuous appearance of a non zero 
average slope below the critical point. In order to charac- 
terize this phase transition quantitatively, we introduce 
an order parameter (a), representing net average slope. 
Results obtained from numerical simulations of the order 
parameter (a), as well as the order parameter suscepti- 
bility and the penetration depth of the slope will be given 
in Section III. 

We also study fractal properties of the avalanches by 
calculating their fractal dimension in the entire range of 
values of the parameter p (Section IV). Our results sug- 
gest that for values of p above the transition point the 
avalanches can be described as self-affine fractals, i.e. the 
shape of the avalanches exhibits an anisotropic behavior. 
Below the transition point both properties, the fractality 
and anisotropy of the relaxation clusters, are lost. 

In Section V we present a detailed analysis of the prob- 
ability distributions of the avalanches and address the 
question whether the model exhibits SOC for p < 1. A 
finite size scaling analysis of the probability distributions 
and a check of certain scaling equations yielded that the 



II. MODEL 

We consider a two-dimensional sandpile model on a 
square lattice of size L x L and integer variables h(i,j) 
representing the local height. We assume a directed 
dynamics, i.e. particles are restricted to flow in the 
downward direction (increasing i). According to the 
widespread 'sandpile language', the first row (i = 1) and 
the last row (i = L) represent the top and the bottom 
of the pile, respectively. To minimize the influence of 
the horizontal boundaries we limit our investigations to 
periodic boundary conditions in this horizontal direction 
(j-direction) . Any site of the lattice has two downward 
and two upward next neighbours, namely h(i+l,j±) and 
h(i - 1, j±) with J± = j ± |(1 ± (-1)*). 

We perturb the system by adding particles at a random 
place on the top of the pile according to 

h(l,j) i — ► /i(l,j') + l, with random j. (1) 

A site is called unstable if the height h(i,j) or at least 
one of the two slopes 

ir(i,3±) = h(i,j)-h(i + l,j ± ) (2) 

exceeds its respective critical value, i.e. if h(i,j) > h c or 
v(hj±) > °c- 

If both slopes are critical one particle after another 
drops alternatively to the next downward neighbours 
until both slopes become subcritical. In the case that 
one slope exceeds a c particles drop to the corresponding 
downward neighbour. 

In contrast to the critical slope, the critical height con- 
dition has a stochastic character. If the local height ex- 
ceeds the critical value h c toppling occurs only with the 
probability p, and then two particles drop to the two 
downward neighbours. 

Because each relaxing cell changes the heights or slopes 
of its four neighbouring sites, the stability conditions are 
applied at these four activated sites in the next updating 
step. Toppling may take place after adding a particle on 
the first row. We apply the slope and the height stabil- 
ity condition simultaneously for all activated sites. An 
avalanche stops if all activated sites are stable. Then we 
start again according to Eq. ([!]). 

It should be emphasized that, in addition to the insta- 
bility criteria, a site is considered as potentially unstable 
only when a particle drops at or near the site. Owing to 
the probabilistic character of the critical height toppling 
rule, locally heights could remain above critical after one 
relaxation event. In this respect, the present model for 
< p < 1 is different from already known critical height 
EIHP and critical slope models. 

One can interpret the parameter 1 — p as being due to 
static friction between the sand grains, which prevents 
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toppling even if the height exceeds the critical value. In 
the case p = 1, our model is identical to the model of 
Dhar and Ramaswamy 0, which exhibits a robust SOC 
behavior. In the limit p — 0, our directed rules lead to 
a steady state where all slopes are equal, i.e. cr(i,j±) — 
<j c — 1. Thus, if a particle is added at the first row it 
performs a directed random walk downward the pile until 
it reaches the boundary and drops out of the system — 
no SOC can occur in this limit. 

We fix the critical height h c and the critical slope <j c 
and restrict ourselves to h c = 2 and a c = 8, although 
it should be emphasized that the results depend on the 
choice of these parameters. 

The behavior of our model is characterized by 
anisotropy effects caused by the preferred direction of 
the assumed dynamics. This anisotropy may influence 
for instance the fractal properties of the avalanches. If 
a measured quantity depends on the direction we intro- 
duce in these cases two different values which describe 
the behavior in parallel and perpendicular direction re- 
spectively, e.g. the fractal dimensions d» and d±. The 
parallel direction corresponds to the preferred direction 
of the dynamics (increasing i). 



III. NONEQUILIBRIUM PHASE TRANSITION 

In this section we consider the general behavior of the 
system, i.e. how the transition from the CHM to the 
CSM takes place. As mentioned above a CHM is char- 
acterized by a zero average slope in contrast to a CSM 
where the steady state displays a non zero average slope. 
In order to describe this transition quantitatively, we in- 
troduce as an order parameter the mean slope of the 
system M 



1 



(3) 



h3 



where 

events, and a. 
defined as 



.) stands for the average over total number of 
n is the local slope at site which is 



(4) 



The summation in Eq. (g) runs from j = 1, 2, 3, L in 
perpendicular direction, and from i = 21, 22..., L — 20 in 
parallel direction, in order to suppress the boundary ef- 
fects from the first and last row of the system. Of course 
the obtained results are independent of the value of the 
cut off length, provided that it is large enough. To nor- 
malize the order parameter one has to choose Ln = L — 40 
and L± = L. The nature of the transition as well as 
the transition point p c should be determined from the 
p-dependence of the average slope Eq. (^) in the vicinity 
of p c - 
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FIG. 1. Average slope of the pile (a) plotted vs. 
probability p. The solid line inside the inset corre- 
sponds to a fit according to Eq. ^. 

In Fig. [j] we show the results of numerical simulations 
of the order parameter defined in Eq. (^). The order 
parameter vanishes at the transition point as, i.e., 



(<r) ~ (Pc ~pY 



(5) 



as one can see from the inset of Fig. [!} We determine the 
transition point p c — 0.293 ±0.002 and the exponent (3 = 
0.8 ± 0.05 for automaton sizes L = 64, 96, 128, 160, 192 . 

We next analyze the fluctuation of the order parame- 
ter. An ordering susceptibility can be defined as 



X = L\\L ± ((a 2 ) -(a) 2 ). 



(6) 



Unfortunately the average slope has as a self-averaging 
character in the critical height regime, that means x van- 
ishes for L — > oo due to a cancellation of the correlations 
(&ij(Tkl) in the interior of the pile. Consequently the sus- 
ceptibility depends only on the boundary term and is of 
relative magnitude ~ L7 . In the case p = 1 the order 
parameter susceptibility can be calculated exactly and is 
given by 



X = 



1 

2Ln 



(7) 



In Fig. |2j we display the susceptibility \ for varying p 
and four different automata sizes. For all considered au- 
tomata sizes x is sharply peaked at the transition point 
p c = 0.293. 

In the inset of Fig. || we plot x ' ^|| vs - P- All different 
curves collapse, except of the values which are close to 
the transition point. The susceptibility scales as 



X~Lf(p-p )-r 



(8) 



with 7 = 0.98 ± 0.05. One can see from Fig. ^ that 
this power law behavior of x holds away from p c up to a 
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certain cut off value which tends to p c if L — > oo. Due 
to this L-dependence the order parameter susceptibility 
tends to zero with increasing L for p > p c . 



of p. For p < 1 a finite slope occurs at the boundaries 
and the pile grows up from the bottom of the pile. The 
slopes penetrate the hole system from the edge of the 
automaton. 
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FIG. 2. The p-dependence of the order parameter 
susceptibility for different system sizes L. The inset 
displays the scaling bahavior of the height regime. 
The solid line corresponds to a power law with expo- 
nent 7 = 0.98. 

Below the critical point our data suggest that the sus- 
ceptibility displays a complicated scaling behavior (see 
Fig. ^). ft is still an open question whether x diverges 
with increasing systems size L or converges to a finite i.e. 
L-independent value. To address this question one has 
to simulate extremely large systems which is beyond our 
present computer capacity. 

In order to break off the self-averaging character of the 
average slopes we defined a nonlinear-susceptibility 



X nl =L\\L ± ((<,>- (<r n! ) 2 ) 
with 



j 2 -. 

L n L± ^ 13 

11 hi 



(9) 



(10) 



We found that x n i shows a cusp at p c , i.e. the nonlinear- 
susceptibility is independent of L in both regimes (not 
shown). This behavior is consistent with the exact value 
for p = 1 : 



The non-diverging behavior of the susceptibilities in- 
dicate that the non equilibrium phase transition is not 
characterized by a diverging correlation length. Thus we 
address the question how the finite slope occurs by ap- 
proaching the transition point. In Fig. ^| we display the 
average height profile (h)(i) of the pile for certain values 




FIG. 3. The average height profile (h(i)) along the 
parallel direction. The penetration depth A is deter- 
mined as the distance from the maximum height to 
the right boundary of the system. The inset displays 
the p-dependence of A. 

We determine the penetration depth A of the slopes 
by measuring the distance of the profile maximum from 
the right boundary (see Fig. 0). The obtained results are 
plotted in the inset of Fig. || for automaton sizes L = 
100, 200. Close to the transition point the penetration 
depth obeys a power law behavior like 

A~(p-p c )-"». (12) 
Our numerical analysis yielded — 0.86 ± 0.06. 



IV. FRACTAL PROPERTIES OF THE 
AVALANCHES 

In this section we consider the scaling properties of the 
avalanches. We define the size s of the avalanches as the 
number of toppled sites in one event. The length of an 
avalanche is defined as the distance from the first row to 
the toppled site with maximal value of the index i. Note 
that this length is by definition parallel to the preferred 
direction, and will be denoted as In. We introduce the 
average width l± of an avalanche via 



(13) 



where l±(i) is the width of an avalanche in the ith row. 
Since there is some arbitrariness in the definition of l± 
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we both tried q = 1 and q = 2 with the result that the 
scaling properties of the width are practically the same in 
both cases. For the sake of simplicity we present in this 
paper only the results for q = 2. We assume that the size 
of avalanches is scale invariant both with the length and 
width. In order to take the anisotropy of avalanches into 
account, we introduce two different scaling exponents d|| 
and d±, corresponding to different directions of scaling in 
and l± , respectively Jl(],[ll| ■ We assume that the average 
size (s) of all avalanches of length hi and of width l±, 
respectively, scales as 



and (s)i ± ~ l±° 



(14) 



If both scaling relations are valid the width scales with 
the length like 



Z x ~^, with c = -A 

01 



(15) 



where the exponent £ describes the anisotropic behavior 
of the avalanches. Since there are only two independent 
exponents we measured d\\ and £ according to Eq. Jl3| ) 
-Eq. @. 

We find that with decreasing p the parallel fractal 
dimension dii increases from the exactly known value 
<i|| = | in the limit p — 1 to two at the transition 
point. In the critical slope regime d\\ remains constant 
and d {l = 2 g. 




L=100 



100 



FIG. 4. Double logarithmic plot of the average 
width vs. length for various values of p as indi- 
cated. All curves obey a power law behavior (solid 
lines) with different accuracy. The inset displays the 
p-dependence of the anisotropic exponent £ obtained 
from the fits according to Eq. [L5| 

Fig. H shows how the width depends on the length for 
different values of p. In all cases we found that Eq. (|l|) 
is valid and extracted values of the anisotropy exponent 



£ (see inset to Fig. ^). With decreasing p the anisotropy 
exponent £ increases from £ = i at p = 1 and reaches 
its maximum value £ ~ 1 for p < p c . We expect that 
this expression is exactly given by £ = 1 and that the 
observed deviations in the critical slope regime are caused 
by finite size effects. The origin of the huge error bars for 
p < p c is a finite curvature of the corresponding curves in 
the log-log plots. These curvatures would be reduced by 
studying larger systems leading eventually to the result 
£ = 1 with higher accuracy. Thus, the critical height 
regime is characterized by an anisotropic shape of the 
avalanches (£ < 1). This behavior changes to an isotropic 
shape (£ = 1) on approaching the transition point. Due 
to the dependence of the scaling factors on the directions 
the avalanches can be regarded as an example of self- 
affine fractals. In analogy to self-similar fractals, it is 
possible to describe self-affine fractals by single exponent 
d f , given by fll[] 



d f = 2- 



£ + 1 



(16) 




0.2 0.4 0.6 

P 



FIG. 5. Fractal dimension df of the avalanches for 
various values of p. 

In Fig. H we show the results for the fractal dimension 
df versus p in the critical height regime. As seen from 
Fig. [|, the relaxation clusters exhibit a fractal behavior 
(df < 2) in the critical height regime p > p c , with ex- 
ception of the case p = 1 . Although the deviation from 
df = 2 is very small they cannot be explained by sta- 
tistical fluctuations due to the size of the error bars and 
the systematic way of the deviations from 2. Below the 
transition point, i.e., for p < p c , the fractal dimension 
df = 2, corresponding to compact clusters. 

One can easily explain this behavior by considering the 
shapes of the avalanches (see Fig. 4 in jsj). In the pure 
critical height model (p = 1) the avalanches are compact, 
i.e. no holes can occur in this limit. With decreasing p 
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some supercritical sites may remain in the interior of the 
avalanche due to the stochastic character of the dynam- 
ics, leading to holes and branching processes. This results 
in a fractal character of the avalanches. By approaching 
the transition point p c this structure of the avalanches is 
lost. In this region the dynamics is more and more char- 
acterized by toppling processes due to the critical slope 
condition. The slope between a toppled site and its two 
backward neighbours may exceed the critical value and 
toppling occurs again. Owing to this multiple topplings 
the probability to find holes inside an avalanche decreases 
with decreasing p and eventually vanishes at p c . At the 
same time branching of the relaxation clusters become 
more rare. 

In this way, we may conclude that the phase transition 
from the critical height to the critical slope regime is ac- 
companied by continuous change from fractal avalanches 
with an anisotropic shape above p c to compact and 
isotropic avalanches below p c . 

V. AVALANCHES DISTRIBUTIONS 

In this section we examine how the SOC behavior of 
the Dhar model vanishes with decreasing p. In the crit- 
ical state, due to the absence of a characteristic length 
scale, the probability P(l\\) for an avalanche longer than 
a certain length lu obeys a power-law, i.e., 

P(l\0~lf K . (17) 

Similarly, the probability distribution of an avalanche size 
larger than s exhibits a power law with another exponent 
via 

P(s) ~ s-\ (18) 

where s is the number of different sites that relaxed in 
one event. Multiply relaxed sites are counted only once 
in this distribution. Double logarithmic plots of the dis- 
tributions P(l\\) and P(s) for various values of the pa- 
rameter p and L = 200 are shown in Fig. || and Fig. 

TABLE I. The exponents for different values of toppling 
parameter p. Due to a lack of a power law behavior of the 
corresponding quantities the exponents r, « and 9 are not de- 
fined for p < 0.5. 
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FIG. 6. Double logarithmic plot of the probabil- 
ity distribution P(l\\) of avalanches of length for 
various values of p as indicated. The solid lines cor- 
respond to power laws. 

One can see from Fig. ^| that the probability distribu- 
tion P(l\\ ) displays no cut-off in the critical height regime. 
For p < p c a sharp cut-off length occurs in the distri- 
bution P (iii). The dynamics is now dominated by the 
toppling processes due to the critical slope condition and 
thus the edge of the system (t = L) influences the top- 
plings at a certain number of rows close to the edge. This 
behavior is typical for critical slope models. 
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FIG. 7. Double logarithmic plot of the probability 
distribution P(s) of avalanches of size s for various 
values of p as indicated. The solid lines correspond 
to power laws. 

Power law behavior according to Eq. |l7| and Eq. [l8| 
occurs only in the critical height regime. The critical 
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exponents k and r are determined by the slopes of the 
linear sections of these curves for p > 0.5. For lower val- 
ues of p we found a finite curvature in the log-log plots of 
the distributions, indicating that the power-law behavior 
is lost for a range of values of p preceding the transi- 
tion point p c . The obtained values of the exponents for 
certain values of p are shown in Fig. || and are listed in 
Table |j| In the limit p = 1 the numerical values of the 
exponents coincide with the exact values [Q k 



i and 



0.90 - 



ID 

g 0.70 
c 

m 0-60 
CO 

o 



O 0.50 - 



' i i 



'4 



• T 

0.40 ♦ K 



I d„ X 



0.30 



0.40 0.50 0.60 0.70 0.80 0.90 1.00 
P 

FIG. 8. Critical exponents « and r determined 
from the power law fits according to Eq. ^ and 
Eq. |l^. In the case of the exponent r the error bars 
are smaller than the symbols. The scaling equation 
Eq. |^ holds for p > 0.6. 

We should emphasize that the above analysis is not a 
sufficient evidence that the model exhibits SOC for p > 
0.5. To improve the accuracy of the results we analyse 
the effects due to the finite size of the systems by using 
simple finite-size scaling M: 



P(s,L) = lr 



sL~ 



(19) 



In order that Eq. ( |l8| ) is recovered for large L the uni- 
versal function g{x) must obey a power-law behavior for 
x —> and the three exponents /3 S , v s and r fulfill the 
relation 



(20) 



The scaling exponent v a can be derived from the fractal 
dimension d\\ : For each value of p we find a cut-off of the 
distribution taking place at a size s max which depends on 
L. If finite-size scaling works all distributions P(s, L), in- 
cluding their cut-offs, have to collapse, i.e. the argument 
of the universal function g has to be constant: 



On the other hand we know from Eq. ( |l4| ) that (s)(^||) 
scales with the length in. Thus (s max ) scales with L in 
the same way and both exponents are identical. In this 
way we have calculated the scaling exponents 



(3 S = d» T anu u s = U| 

from the measured exponents d\\ and r. 



and 



(22) 
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FIG. 9. Finite size scaling analysis of the avalanche 
size distribution P(s) for p = 1.0,0.9,0.8,0.7,0.6,0.5 
and five different system sizes L. The curves cor- 
responding to p < 1 are shifted in the upper right 
direction. The finite size scaling ansatz works quite 
well for p > 0.6. 

In Fig. H we show the results of the scaling plot for 
p > 0.5 and five different automata sizes. Finite-size 
scaling works well for p > 0.7. In the case of p < 0.6 
the deviations grow, i.e., the different curves do not com- 
pletely collapse. 

In the following we will proof a scaling relation among 
the exponents r, n and d|| , indicating that the exponents 
of the distributions P(s) and P(l\\) are not independent. 
The relation 



p(h) d h - p(h) 



dZ || 
d7 



ds 



p(s) ds 



and Eq. (n_J) lead to the scaling relation 



(23) 



(24) 



Smax L 8 — COllSt. 



(21) 



Here p(l\\) and p{s) represent the probability density of 
an avalanche of length l\\ and size s, i.e. they obey the 
power-laws p (h \ ) ~ lu~~ K_ andp(s) ~ s _T_1 , with k and 
r defined in ( pLT|) and (|l8|), respectively. 

In Fig. |^ the product d|| r is also shown and compared 
with the values of the exponent n for various values of 
p. We conclude that the scaling relation (E4J) holds for 
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p > 0.7. These results confirm the finite-size scaling anal- 
ysis and suggest that criticality is lost at a certain point 
between p = 0.6 and p — 0.7. We expect that preceding 
the transition point p c the occurrence of toppling pro- 
cesses due to the critical slope condition increases and 
results in a loss of criticality which is connected to the 
pure critical height behavior. 

We have also investigated the avalanches durations 
t and found that the average duration scales with the 
length hi according to 



(t)l 



lu 



(25) 



Analogous to the size s the probability distribution of the 
duration obeys a power-law behavior 



p{t) ~ r 



(26) 



for p > 0.7. In the same area of p finite-size scaling 
works quite well and the corresponding scaling relation 
k = zO is fulfilled. The obtained values of the exponents 6 
and z are listed in Table |. These measurements confirm 
our previous conclusion that the system exhibits SOC for 
p > 0.7 and not in the whole critical height regime. 



VI. CONCLUSIONS 

We have studied numerically a sandpile model with 
stochastic dynamics which exhibits a noncquilibrium 
phase transition as the parameter p (representing static 
friction) is varied. The probability parameter allows us 
to tune the dynamics from a critical height to a critical 
slope behavior. The average net slope (a) plays the role 
of the order parameter and obeys a power law depen- 
dence in p c — p. In contrast to the diverging penetration 
depth both the linear and nonlinear susceptibility shows 
no singular behavior by approaching the transition point 
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FIG. 10. Phase diagram of the model which illus- 
trates the major properties of the investigated system. 
The order parameter a and the penetration depth A 
depend on e = \p — p c \. 



We found that the nonequilibrium phase transition 
is accompanied by a transition from an anisotropic to 
an isotropic behavior of the avalanches. Fractality oc- 
curs only in the case of the anisotropic shape of the 
avalanches. 

Finite size scaling analysis works for the probability 
distributions of the avalanches and scaling relations hold 
for p > 0.7 indicating that the system exhibits SOC in 
this region. The corresponding exponents have a nonuni- 
versal i.e. p-dependent behavior. 

Due to long computation time, we did not simulate lat- 
tices larger than L = 500 in order to check if the power- 
law behavior is still maintained for values of p < 0.6. The 
problem of disappearance of the SOC on approaching the 
nonequilibrium phase transition remains for future study. 

A complete illustration of the behavior of our sandpile 
model is depicted in Fig. [ic| . 
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